Risk and Resilience in Deep Time

Analysis Scripts for Paper Published in Holocene

Set up libraries

Load data

c14dates <- read_csv("c14dates.csv")

Rows: 1452 Columns: 16── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (10): industry, site, unit, calib.curve, sample.id, material, delta 13, method&id, comments, citation
dbl  (6): ID, Use, BP.cal.median, C14.mean, C14.SD, C14.COV
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Figure 1: Journal Publications

rr.papers %>% 
  ggplot(aes(x=years)) + 
  geom_path(aes(y=JAS.risk/JAS.papers, color='Journal of Archaeological Science', lty='risk'), group=1, lwd=1.5) + 
  geom_path(aes(y=JAS.resilience/JAS.papers, color='Journal of Archaeological Science', lty='resilience'), group=1, lwd=1.5) + 
  geom_path(aes(y=Antiquity.risk/Antiquity.papers, color='Antiquity', lty='risk'), group=1, lwd=1.5) + 
  geom_path(aes(y=Antiquity.resilience/Antiquity.papers, color='Antiquity', lty='resilience'), group=1, lwd=1.5) + 
  scale_y_continuous(labels = scales::percent) + 
  scale_color_manual(values = c('Journal of Archaeological Science' = 'red', 'Antiquity' = 'blue')) + 
  scale_linetype_manual(values = c('risk' = 'solid', 'resilience' = 'dashed')) + 
  labs(title = "Papers Mentioning Risk or Resilience", 
       subtitle = "Change in Past 40 Years", 
       x="5-year intervals", 
       y="% of papers published", 
       color='journal:', 
       linetype='topic:') + 
  theme_bw(base_size = 30) + 
  theme(legend.position = 'bottom', legend.key.width =  unit(50, "points"))

Modeling Biocultural Interaction

Figure 2: NGRIP2 and GISP2 entire cores

icecores %>%
  ggplot() + 
  geom_point(aes(x=years.BP, y=d18O.NGRIP2.ppt), color='blue', size=.5, alpha=.2) + 
  geom_point(aes(x=years.BP, y=d18O.GISP2.ppt), color='pink2', size=.5, alpha=.2) + 
  geom_smooth(aes(x=years.BP, y=d18O.NGRIP2.ppt, color='ngrip2', lty='ngrip2'), method='loess', span=.01, se=FALSE, lwd = 1)  + 
  geom_smooth(aes(x=years.BP, y=d18O.GISP2.ppt, color='gisp2', lty='gisp2'), method='loess', span=.01, se=FALSE, lwd = 1) + 
  scale_colour_manual(name=NULL, values =c('gisp2'='indianred3','ngrip2'='navy'), labels = c('GISP2','NGRIP2')) + 
  scale_linetype_manual(name=NULL, values =c('gisp2'='solid','ngrip2'='solid'), labels = c('GISP2','NGRIP2')) + 
  scale_x_reverse() + 
  #scale_x_continuous(breaks=c(19000,14000,10000), labels = c("19000","14000","10000"), trans = "reverse") + 
  geom_vline(xintercept = c(27000,22000,11000), lty='dashed', lwd=1) + 
  labs(x="years cal BP", 
       y="delta O18\ncolder << — >> warmer", 
       color="Cores",
       title="Paleoclimate Proxies from Greenland Ice Cores") +
  theme_minimal(base_size = 18) + 
  theme(legend.position="bottom", 
        plot.title = element_text(face = "bold", hjust = 0.5), 
        )

Figure 3: Box Plots with Condensed Fitness Categories

Hominin_demography %>% 
  mutate(fitness.aggregated = Fitness, 
         fitness.aggregated = recode(fitness.aggregated, 
                                     "birthrate=deathrate"="equal fitness", 
                                     "birthrate>deathrate"="equal fitness", 
                                     "MM higher mortality"="MM less fit", 
                                     "MM lower fertility"="MM less fit", 
                                     "MM lower mortality"="MM more fit", 
                                     "MN higher mortality"="MN less fit", 
                                     "MN lower fertility"="MN less fit", 
                                     "MN lower mortality"="MN more fit", 
                                     "NN higher mortality"="NN less fit", 
                                     "NN lower fertility"="NN less fit", 
                                     "NN lower mortality"="NN more fit"),
         fitness.aggregated = factor(fitness.aggregated, 
                                     levels = c(
                                       "equal fitness", 
                                       "MM less fit", 
                                       "MM more fit", 
                                       "MN less fit", 
                                       "MN more fit", 
                                       "NN less fit", 
                                       "NN more fit"))) %>% 
  ggplot(aes(x=phenotype, y=population)) + 
  geom_boxplot(fill="grey") + 
  facet_grid(fitness.aggregated~factor(home.range.radius), scales="free_y") + 
  labs(x="agent phenotype", y="population at end of simulation") + 
  theme_bw(base_size = 16)

Figure 4: Example Model Trajectory with Logistical Foraging (Foraging Radius=12)

Hominin_demography_trajectory %>% 
  ggplot(aes(y=population, x=step)) + 
  geom_line(aes(color=phenotype), size=2) +
  scale_color_manual(values = c("red", "orange", "grey", "forestgreen", "blue")) + 
  labs(x="model step", y="population", color="agent phenotype") + 
  theme_bw(base_size = 16) + 
  theme(legend.position="bottom")

Demographic Dynamics of West Mediterranean

Figure 5a: Summed Probability Distribution (SPD) of Radiocarbon Dates in the West Mediterranean

Calibrate dates

c14dates$BP.cal.median <- c14dates %>% 
  with(., calibrate(x=C14.mean, errors=C14.SD, calCurves = calib.curve, normalised=TRUE, calMatrix=FALSE)) %>% 
  medCal()

Calculate SPD for entire region

Dates file filtered to remove all date with COV > 0.1

Graph SPD

par(mar=c(7,10,7,3))
plot(model.all, xlim = c(35000,8000), col.obs = 'red', lwd.obs = 3, drawaxes = F)
axis(1, cex.axis = 1, pos = -.0006)

axis(2, cex.axis = 1, pos = 35500)
mtext(side=1, line=3, "years cal BP", cex=1.3)

mtext(side=2, line=4, "summed probability", cex=1.3)
title(main=paste("West Mediterranean Late Pleistocene Assemblages (N = ", nrow(subset(c14dates, calib.curve != "normal")), ")"), cex.main = 1.5)

Figure 5b: SPD Permutation Test of Aggregated Technocomplexes

Classify aggregated technocomplexes

c14dates <- c14dates %>% 
  mutate(industry.group = 
           case_when(
             str_detect(tolower(industry), "aurig") ~ "Aurignacian", 
             str_detect(tolower(industry), "gdrav") ~ "Gravettian",  
             str_detect(tolower(industry), "solut") ~ "Solutrean",
             str_detect(tolower(industry), "middle solutrean") ~ "Solutrean", 
             str_detect(tolower(industry), "salp") ~ "Solutrean",              
             str_detect(tolower(industry), "epig") ~ "Magdalenian",  
             str_detect(tolower(industry), "epi-gr") ~ "Magdalenian",
             str_detect(tolower(industry), "eppig") ~ "Magdalenian",
             str_detect(tolower(industry), "eppgr") ~ "Magdalenian",
             str_detect(tolower(industry), "apig") ~ "Magdalenian",
             str_detect(tolower(industry), "magd") ~ "Magdalenian",
             str_detect(tolower(industry), "badeg") ~ "Magdalenian",
             str_detect(tolower(industry), "gravet") ~ "Gravettian",  
             str_detect(tolower(industry), "gravt") ~ "Gravettian",  
             str_detect(tolower(industry), "azil") ~ "Epipaleolithic", 
             str_detect(tolower(industry), "romanell") ~ "Epipaleolithic", 
             str_detect(tolower(industry), "epipal") ~ "Epipaleolithic", 
             str_detect(tolower(industry), "mesol") ~ "Mesolithic", 
             str_detect(tolower(industry), "sauve") ~ "Mesolithic", 
             str_detect(tolower(industry), "tarden") ~ "Mesolithic", 
             str_detect(tolower(industry), "montcl") ~ "Mesolithic", 
             str_detect(tolower(industry), "montad") ~ "Mesolithic", 
             str_detect(tolower(industry), "castel") ~ "Mesolithic"
             ), 
         industry.group = factor(industry.group, 
                                 levels = c("Aurignacian", 
                                            "Gravettian", 
                                            "Solutrean", 
                                            "Magdalenian", 
                                            "Epipaleolithic", 
                                            "Mesolithic")),
         .after=industry)

Permutation test of aggregated technocomplexes

Permutation test results

summary(perm.test)

'permTest()' function summary:

Number of sets: 6
Statistical Significance computed using 200 simulations. 

--- Aurignacian ---
Number of radiocarbon dates:94
Number of bins:84

Global p-value: 0.00498

Significant positive local deviations at:
40000~29265 BP 
28924~28146 BP 
27436~26419 BP 

Significant negative local deviations at:
18435~12504 BP 
12048~11097 BP 
10230~7646 BP 

--- Gravettian ---
Number of radiocarbon dates:131
Number of bins:109

Global p-value: 0.00498

Significant positive local deviations at:
32846~25171 BP 
24609~24371 BP 
24318~24295 BP 
24181~24163 BP 

Significant negative local deviations at:
21508~21288 BP 
20184~19623 BP 
16356~13859 BP 
13809~10550 BP 
9793~7266 BP 

--- Solutrean ---
Number of radiocarbon dates:120
Number of bins:88

Global p-value: 0.00498

Significant positive local deviations at:
26079~19630 BP 

Significant negative local deviations at:
31656~31025 BP 
18570~18569 BP 
18560~18229 BP 
18121~18117 BP 
13728~12489 BP 
12168~7419 BP 

--- Magdalenian ---
Number of radiocarbon dates:516
Number of bins:350

Global p-value: 0.00498

Significant positive local deviations at:
20402~20034 BP 
19782~12264 BP 

Significant negative local deviations at:
36464~24702 BP 
24446~23560 BP 
23394~23225 BP 
23221~22834 BP 
11441~10975 BP 
10860~7250 BP 

--- Epipaleolithic ---
Number of radiocarbon dates:120
Number of bins:94

Global p-value: 0.00498

Significant positive local deviations at:
13305~10489 BP 
10483~10341 BP 
10146~9669 BP 

Significant negative local deviations at:
32202~26954 BP 
26507~25629 BP 
25500~22554 BP 
20524~19706 BP 
18871~17393 BP 
16130~14915 BP 

--- Mesolithic ---
Number of radiocarbon dates:312
Number of bins:206

Global p-value: 0.00498

Significant positive local deviations at:
11324~7000 BP 

Significant negative local deviations at:
34855~11918 BP 

Graph Permutations

plot(perm.test, 
     focalm = "Aurignacian", 
     xlim = c(35000,8000), 
     col.obs = 'red', 
     lwd.obs = 3,
     bbty = 'b', 
     ylim = c(0, 0.04),
     drawaxes = F)
  axis(1, cex.axis = 1, pos = -.0006)

  axis(2, cex.axis = 1, pos = 35500)
  mtext(side=1, line=3, "years cal BP", cex=1.3)

  mtext(side=2, line=4, "summed probability", cex=1.3)
  title(main=paste("Aurignacian (N = ", nrow(subset(c14dates, industry.group == "Aurignacian")), ")"), cex.main = 1.5)

 
plot(perm.test, 
     focalm = "Gravettian", 
     xlim = c(35000,8000), 
     col.obs = 'red', 
     lwd.obs = 3,
     bbty = 'b', 
     ylim = c(0, 0.04),
     drawaxes = F)
  axis(1, cex.axis = 1, pos = -.0006)

  axis(2, cex.axis = 1, pos = 35500)
  mtext(side=1, line=3, "years cal BP", cex=1.3)

  mtext(side=2, line=4, "summed probability", cex=1.3)
  title(main=paste("Gravettian (N = ", nrow(subset(c14dates, industry.group == "Gravettian")), ")"), cex.main = 1.5)

  
plot(perm.test, 
     focalm = "Solutrean", 
     xlim = c(35000,8000), 
     col.obs = 'red', 
     lwd.obs = 3,
     bbty = 'b', 
     ylim = c(0, 0.04),
     drawaxes = F)
  axis(1, cex.axis = 1, pos = -.0006)

  axis(2, cex.axis = 1, pos = 35500)
  mtext(side=1, line=3, "years cal BP", cex=1.3)

  mtext(side=2, line=4, "summed probability", cex=1.3)
  title(main=paste("Solutrean (N = ", nrow(subset(c14dates, industry.group == "Solutrean")), ")"), cex.main = 1.5)

  
plot(perm.test, 
     focalm = "Magdalenian", 
     xlim = c(35000,8000), 
     col.obs = 'red', 
     lwd.obs = 3,
     bbty = 'b', 
     ylim = c(0, 0.1),
     drawaxes = F)
  axis(1, cex.axis = 1, pos = -.0006)

  axis(2, cex.axis = 1, pos = 35500)
  mtext(side=1, line=3, "years cal BP", cex=1.3)

  mtext(side=2, line=4, "summed probability", cex=1.3)
  title(main=paste("Magdalenian (N = ", nrow(subset(c14dates, industry.group == "Magdalenian")), ")"), cex.main = 1.5)

  
plot(perm.test, 
     focalm = "Epipaleolithic", 
     xlim = c(35000,8000), 
     col.obs = 'red', 
     lwd.obs = 3,
     bbty = 'b', 
     ylim = c(0, 0.04),
     drawaxes = F)
  axis(1, cex.axis = 1, pos = -.0006)

  axis(2, cex.axis = 1, pos = 35500)
  mtext(side=1, line=3, "years cal BP", cex=1.3)

  mtext(side=2, line=4, "summed probability", cex=1.3)
  title(main=paste("Epipaleolithic (N = ", nrow(subset(c14dates, industry.group == "Epipaleolithic")), ")"), cex.main = 1.5)

  
plot(perm.test, 
     focalm = "Mesolithic", 
     xlim = c(35000,8000), 
     col.obs = 'red', 
     lwd.obs = 3,
     bbty = 'b', 
     ylim = c(0, 0.12), 
     drawaxes = F)
  axis(1, cex.axis = 1, pos = -.0006)

  axis(2, cex.axis = 1, pos = 35500)
  mtext(side=1, line=3, "years cal BP", cex=1.3)

  mtext(side=2, line=4, "summed probability", cex=1.3)
  title(main=paste("Mesolithic (N = ", nrow(subset(c14dates, industry.group == "Mesolithic")), ")"), cex.main = 1.5)